Statistical modeling of trends in infant mortality after atmospheric nuclear weapons testing

The global fallout from atmospheric nuclear weapons testing in the 1950s and 1960s caused by far the greatest exposure of mankind to ionizing radiation. Surprisingly few epidemiological studies of the possible health effects of atmospheric testing have been conducted. Here, long-term trends in infant mortality rates in the United States (U.S.) and five major European countries (EU5) were examined: The United Kingdom, Germany, France, Italy, and Spain. Bell-shaped deviations from a uniformly decreasing secular trend were found beginning in 1950, with maxima around 1965 in the U.S. and 1970 in EU5. From the difference between observed and predicted infant mortality rates, in the period 1950–2000, the overall increase in infant mortality rates was estimated to be 20.6 (90% CI: 18.6 to 22.9) percent in the U.S. and 14.2 (90% CI: 11.7 to 18.3) percent in EU5 which translates to 568,624 (90% CI: 522,359 to 619,705) excess infant deaths in the U.S. and 559,370 (90% CI: 469,308 to 694,589) in the combined five European countries. The results should be interpreted with caution because they rely on the assumption of a uniformly decreasing secular trend if there had been no nuclear tests, but this cannot be verified. It is concluded that atmospheric nuclear weapons testing may be responsible for the deaths of several million babies in the Northern Hemisphere.


Introduction
The testing of nuclear weapons in the atmosphere resulted in worldwide exposure to radioactive fallout from the explosions. Atmospheric nuclear testing was at its height in the late 1950s and early 1960s. The cumulative explosive power of the tests corresponded to 545 megatons of 2,4,6-trinitrotoluene (TNT), which is equivalent to 40,000 atomic bombs of the type that was dropped on Hiroshima at the end of the Second World War (WW2). The collective dose to the world population was estimated at 30 million Person-Sievert (PersSv) [1] p.99), which can be compared with 600,000 PersSv from the Chernobyl accident in 1986 [1] p.115). Fig 1 shows the annual strontium deposition in the northern hemisphere from atmospheric nuclear testing. Internal exposure was greatest in the Northern Hemisphere where most of the weapons testing took place.
Epidemiological studies of possible health effects from atmospheric testing focused on childhood leukemia. Radiation-induced excess risk of childhood leukemia is apparent within a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 some five years of exposure, so the increase in childhood leukemia after exposure to fallout is expected to appear in the late 1960s. A temporal correlation study of childhood leukemia in some Nordic countries (Denmark, Finland, Iceland, Norway, and Sweden) with the fallout from atmospheric nuclear weapons testing found little evidence of increased incidence of leukemia among children born in these years [3]. During the high exposure period, children would have been subjected to an additional equivalent dose of around 1.5 mSv, similar to doses received by children in several parts of central and eastern Europe owing to the Chernobyl accident, and about 50% greater than the annual dose equivalent to the red bone marrow of a child from natural radiation [3]. Rates of leukemia in the high exposure period were slightly higher than in the surrounding medium exposure period (relative risk for ages 0-14: 1.07, 95% confidence interval 1.00 to 1.14; for ages 0-4: 1.11, 1.00 to 1.24). A later study also found "no evidence of a wave of excess cases (of childhood leukemia incidence) corresponding to the peak of radioactive fallout from atmospheric weapons testing which was expected based on conventional risk estimates" [4].
In the late 1960s, Sternglass investigated trends in infant mortality rates in the United States following the Nevada test series and the Pacific tests in the 1950s [5]. After 1950, infant mortality rates deviated from an exponentially declining trend observed in 1935-1950 (see S1 Fig in S1 File). Citing early warnings by Sakharov [6], Sternglass attributed the excess infant mortality in the United States to immune system damage from strontium-90 in the fallout from atmospheric nuclear weapons testing. He estimated that in the U.S. "some 400,000 infants of less than one year of age probably had died as the result of nuclear fallout between 1950 and 1965" [7]. Sternglass' claims have been questioned [8][9][10][11]; a major criticism was that "a careful evaluation of the available data does not demonstrate a direct relationship between fallout and the observed trend changes" [10].
Whyte analyzed data on first-day neonatal mortality from England and Wales and the U.S. [12]. The trends were very similar in the two countries: an upward deviation from the trend of the preceding years starting in 1950, with a broad maximum around 1965 and a decline until  about 1980 when the data resumed the trend of the years 1935-1950, see S2 Fig in S1 File. For  the United States, Whyte identified 280,000 excess neonatal deaths during 1951-1980 compared with the trend during the remaining years of the study period ( . Among environmental factors which might explain the excess mortality, Whyte mentions exposure to strontium-90 resulting from atmospheric weapons testing.
The present author investigated the development of perinatal mortality in West Germany in the period 1955-93 and found similar deviations from a monotonously decreasing trend, but with a maximum around 1970 [13], see S3 Fig in S1 File. As a possible explanation for the time lag between exposure and effect, he suggested that radioactive strontium 90 ingested by girls at about age 14, during the period of increased bone growth, impairs their immune systems, which in turn may lead to increased mortality of children born to these girls at later ages.
The present exploratory study adopts the approach used in [13] to examine infant mortality trends in the United States before and after atmospheric nuclear weapons testing and compares the results with those of the combined data from five European countries (United Kingdom, Germany, France, Italy, and Spain).

Data and methods
As mentioned above, Körblein reported an association between perinatal mortality in West Germany and the calculated strontium burden of pregnant women [13]. His rather complex analysis included data on annual strontium concentration in the fallout and the maternal age distributions for each year of the study period . However, a nearly perfect fit to this data was also obtained with a regression model using two bell-shaped excess terms (lognormal density distributions), superimposed on a uniformly decreasing trend, see S4 Fig in S1 File. Therefore, lognormal distributions are used in the present study to model the trend of infant mortality rates. Infant mortality was chosen as the endpoint because such data are available online for many countries, which is not the case for perinatal mortality.
Infant death is defined as the death of a live-born infant within the first year. Infant mortality is the number of infant deaths per 1000 live births. Annual numbers of live births and infant deaths for the United States (U.S.) and the United Kingdom (UK), France, Italy, and Spain are provided at https://mortality.org. German infant mortality data for 1931-1943 were obtained from the German Federal Office of Statistics on request. The German data from 1946 to 2021 are published in an article available for download on their website, https://www.destatis.de [14]. Data for 1944 and 1945 do not exist.
The present analysis of infant mortality in the United States uses data from 1934 through 2018. For Europe, the pooled data from the United Kingdom (UK), Germany, France, Italy, and Spain, 1931-2018, were chosen.
Two competing mathematical models were used for the undisturbed trend: (1) a modified logistic function allowing for a lower limit of infant mortality, as in [13]; and (2) an exponential trend with a third-degree polynomial of time. The full regression Model (1) with 3 bellshaped excess terms is non-linear and has the following analytical form: Here, E(y(t)) is the predicted value of infant mortality y(t). Time t is the calendar year minus 1930; the constant α is the lower limit of infant mortality; β 1 and β 2 are trend parameters; and β 5 through β 13 are parameters defining the three lognormal functions.
Model (2) has the following form: The variables t2 and t3 are time t to the second and third power; β 3 and β 4 are additional trend parameters.
For reweighting, binomial variances var = fit�(1−fit)/LB were used, where fit is the fitted value and LB is the number of live births. Two-sided statistical tests were applied, and p-values <0.05 were considered statistically significant. R statistical software version 4.2.2 was used for regression analyses and plotting graphs [15].
To determine the confidence limits for the estimated number of excess deaths in 1950-2000, the result for the regression parameter β 2 in Model (1) was used. Because the research question is one-sided, i.e., whether there is increase-not decrease-in infant mortality, the 90% confidence interval (90% CI) for the number of excess deaths was determined. The 90% CI of the parameter β 2 in Model (1) was calculated with the R function confint() which uses Monte Carlo simulation to estimate the upper and lower confidence limits. The results of regressions with parameter β 2 fixed to the upper / lower limits were used to calculate the predicted infant deaths and the number of excess cases. The points in time of the maxima of the bell-shaped excess terms were calculated as exp (μ−σ 2 ) where μ and σ are the parameters defining the lognormal function.

Results
To start with, former studies of infant mortality in the U.S. shall be revisited. Sternglass compared the observed infant mortality rates after 1950 with the rates predicted from the extrapolated trend of the data in 1935-1950 to determine any excess infant mortality after the beginning of atmospheric nuclear weapon tests in Nevada in the early 1950s [5]. S5 Fig in S1 File shows the trend in infant mortality in the United States over the period 1934-2018 and the result of regressing the 1934-1950 data on an exponential model. The dashed line is the predicted trend for the period after 1950. For clarity, a semilogarithmic plot is used. Until the end of the study period, the observed infant mortality rates are well above the predicted trend. A purely exponential trend model is therefore inappropriate; it implies that infant mortality moves toward zero with time, which is not the case. As a result, Sternglass' method overestimates the number of excess infant deaths.
Whyte determined the excess mortality in 1951-1980 by interpolation based on data before 1951 and after 1980 using an exponential trend model [12].

United States
The data from the Unites States (U.S.) were first analyzed with Model (1). Regression without bell-shaped excess terms gave deviance = 43,881 (df = 82). With one excess term, the deviance was 3,244 (df = 79); a second excess term reduced the deviance to 2,363 (df = 76). A slight improvement in fit, although not statistically significant (p = 0.205), was achieved when the first excess term was replaced by a superposition of two excess terms with equal half-widths. The final Model (1) fitted the data well, albeit with large overdispersion; the deviance was 2,264 (df = 74), much greater than expected from a random distribution.  Model (2) was then applied to fit the U.S. data. Regression without excess terms gave deviance = 42,396 (df = 81); with one excess term the deviance was 2,998 (df = 78); the addition of a second excess term yielded deviance = 2,414 (df = 75), and the final Model (2) with 3 excess terms reduced the deviance to 2,352 (df = 73). 530,348 excess infant deaths were estimated during 1950-2000, an overall 19% increase. The maxima of the bell-shaped excess terms were found near 1960, 1970, and 1988. Table 1 contains the regression results (parameter estimates with standard errors and t-values) for the two regression models.

PLOS ONE
Infant mortality after atmospheric nuclear weapons testing Table 2 shows the improvement in fit resulting from stepwise increasing the number of excess terms for both regression models. Here, OD (overdispersion) is the deviance divided by the number of degrees of freedom (df2), and df1 is the number of parameters added with each step of model refinement.

Europe
The regression of infant mortality rates in Europe is complicated by the impact of the Second World War (WW2) on infant mortality. Here, the combined data from the United Kingdom (UK), Germany, France, Italy, and Spain, hereafter denoted as EU5, was analyzed from 1950 onward only, assuming that the effects of WW2 on infant mortality had been overcome by 1950. Regression with Model (1) without excess terms yielded deviance = 14,950 (df = 66); with one excess term, the deviance decreased to 2,245 (df = 63); a second excess term with the same width reduced the deviance to 1,642 (df = 61), and regression with the full Model (1) Table 3. Table 4 shows the improvement in fit with stepwise model refinement.  Table 7. The excess infant deaths add up to 547,352.
To estimate any difference in effect size between central and southern European countries, infant mortality data from three central-(UK, Germany, and France) and two southern (Italy and Spain) European countries were analyzed separately with Model (1). The regression results are listed in Table 8. Interestingly, the estimates for the slope (parameter β 2 ) were almost identical, and the estimates of the median values of the main two lognormal distributions agreed within the error bounds (see Table 8, parameters β 6 and β 9 ). In southern Europe, excess infant mortality rates were about twice as high as in central Europe (see Fig 6,

Sensitivity analyses
United Kingdom. The data from the United Kingdom show little overdispersion and lend themselves to more detailed exploratory analysis. Model (1) is supplemented by two additional bell-shaped excess terms. The fourth term reduces the deviance from 259 (df = 58) to 201.6 (df = 55), p = 0.003, and a fifth excess term yields deviance = 143.5 (df = 52), p = 0.0005. Table 9 shows the improvement in fit by gradually increasing the number of excess terms.    The Swedish data was analyzed with the same regression model used above for the UK data, i.e. with Model (1) and five excess terms. The regression results are presented in S1 Table in S1 File, together with the results for the United Kingdom. S9 Fig in S1 File shows infant mortality rates in Sweden together with those from EU5. The mean excess rate in Sweden is smaller than that in EU5 which is consistent with the latitudinal dependence of strontium-90 deposition from global fallout. As in the data from the UK, there is a peak after the Chernobyl accident which may be related to the consumption of contaminated meat from reindeer. The number of excess infant deaths attributable to this spike is 276 in 1986-1994, with a maximum of 127 excess cases in 1991 alone.

Analysis of European data, 1931-2018
To avoid trend bias due to the impact of WW2, the above analyses of European infant mortality data were restricted to data after 1949. As an alternative, the rise in infant mortality during and after WW2 is approximated by an additional bell-shaped excess term. This new model, i.e. Model (2), supplemented by an additional bell-shaped term, will hereafter be referred to as Model (3) Regressions with Model (3) were also conducted for each of the five European countries, see S10-14 Figs in S1 File. To analyze the German data, a second bell-shaped excess term was needed to fit the data before 1950. In France, highly significant peaks in infant mortality were found in 1940 and 1945; these two years were therefore omitted from the regression. The individual regression results for the five European countries are listed in S3, S4 Tables in S1 File. Table 10 shows the number of excess infant deaths in individual countries.
Regression of the Swedish data on infant mortality, 1931-2018, was conducted with Model (3) as well. Highly significant drops in infant mortality were detected in 1942 and 1943 with deviations of more than four standard deviations from the predicted trend; infant mortality

Effect of gender
To check whether the excess infant mortality depends on gender, the data from the U.S.  and from the UK (1950-2018) were analyzed for male and female infants separately with Model (1) and three bell-shaped excess terms. For males and females, the points in time of the maxima (median values) as well as the widths (geometric standard deviations) and the relative excess mortalities agreed within the limits of error for the two main excess terms with peaks around 1960 and 1970. However, the peak in the late 1980s was greater for males than for females. The results for the United States and the United Kingdom are listed in S5 and S6 Tables in S1 File. The respective graphs are shown in S16, S17 Figs in S1 File.

Trend analysis of German stillbirth rates
Since annual data of livebirths and stillbirths from Germany are also provided in [14], a trend analysis of stillbirth rates in the period 1950-2020 was performed using regression Model (2). In 1994, the definition of stillbirth was changed from a birth weight of >1000 grams to >500 grams, so stillbirth rates after 1993 were adjusted by adding two dummy variables that allowed for an increase in 1994 and a level shift in the period 1995-2020. The regression model used a

Infant mortality in Russia
Out of curiosity, data on infant mortality in Russia was also examined. On the mortality.org website, data on live births and infant deaths from Russia are only available from 1959 onwards, so that a meaningful study of the possible effects of nuclear weapons testing is not possible. However, infant mortality rates show a relative maximum in 1976, with a decreasing trend thereafter. An analysis of infant mortality rates from 1976 to 2010 exhibits an upward shift in 1984 and a long-term increase and decrease after the Chernobyl accident in 1986 relative to rates predicted from the pre-Chernobyl trend. The deviation of observed rates from the predicted trend can be approximated by two lognormal distributions with maxima around 1994 and 2001. The excess of infant mortality rates corresponds to 87,436 additional infant deaths during 1987-2010 (see S18 Fig in S1 File).

Discussion
The present study examines the trends in infant mortality in the United States (U.S.), and five major European countries combined (EU5) before and after atmospheric nuclear weapons testing. in the 1950s and early 1960s. In both regions, deviations from a uniformly decreasing trend are found in the period 1950-2000. These deviations were approximated by three superimposed bell-shaped excess terms (lognormal distributions). The median values of the two main lognormal distributions agree within the error limits for the U.S. and the EU5 (see parameters β 6 , β 9 in Tables 1 and 3). This also holds for the comparison of the results for central Europe (UK, Germany, France) and southern Europe (Italy and Spain), see Table 8. The point in time of the peak in excess infant mortality depends on the relative size of the two main excess terms with maxima in the early1960s and 1970s. In the U.S., the first excess term is greater than the second, so the overall maximum appears earlier than in EU5 where the second excess term dominates (cf . Figs 2 and 4).
From the regression results for individual European countries in Tables 5 and 6, the timing of the maxima for the first and second lognormal distributions was calculated, see Table 11. In the United Kingdom, both peaks are lagged compared to the corresponding peaks in EU5, while in Spain, both peaks occur earlier than in EU5. Differences in the age distribution of mothers in the United Kingdom and Spain could explain this observation; Spanish women in the 1960s and 1970s may have been younger at their first pregnancy than women from the United Kingdom, but this could not be verified because the author did not have historical data on the age distribution of mothers.
The two main peaks in infant mortality follow, with a lag of about seven years, the fallout peaks after the two largest atmospheric weapons tests: the U.S. hydrogen bomb "Castle Bravo" at Bikini Atoll, March 1, 1954, with an explosive force equivalent to 15 million tons (Mt) of TNT, and the largest Soviet hydrogen bomb "Tsar Bomba" at Novaya Zemlya, October 30, 1961, with 50 Mt TNT. The dominance of the first peak in the U.S. may be explained by the Nevada tests and higher exposure to fallout from the Pacific tests in the 1950s. Both Sternglass and Whyte had suggested strontium in the fallout as a likely cause of the observed increased infant mortality [5,12,17]. However, no direct relationship between annual strontium exposure and increased infant mortality has been demonstrated [11]. Körblein proposed a possible mechanism to explain the observed long-term deviations of perinatal mortality from a uniform declining trend [13,18]. This mechanism, a delayed effect of strontium on the immune system, was used to model the trend of infant mortality rates in the present study. https://doi.org/10.1371/journal.pone.0284482.g011 Table 11. Estimated timing (calendar year) of the two main infant mortality peaks.

First peak Second peak
Country Peak position (95% CI) Peak position (95% CI) The main limitation of the study is that we do not know how infant mortality trends would have evolved in the absence of the effects of atmospheric weapons testing. Improvements in health care, such as the introduction of prenatal ultrasonography in the early 1980s, and socioeconomic influences likely affect infant mortality trends. The results must therefore be interpreted with due caution. However, the rates of stillbirths in Germany (1950-2020) follow a remarkably smooth trend, with the exception of a bell-shaped increase between 1960 and 1980 (see Fig 10).

Conclusion
The present analysis of infant mortality after atmospheric nuclear weapons testing shows similar trends in the United States and Europe. The same statistical model was used to fit the data in both regions: a uniformly decreasing trend with three superimposed bell-shaped excess terms. The increased infant mortality is interpreted as a delayed effect of strontium on the immune system. However, this hypothesis cannot be tested because there are no data on strontium exposure in pregnant women. Acceptance of the strontium hypothesis means that atmospheric weapons testing could be responsible for the deaths of millions of infants in the Northern Hemisphere.